Regresia segmentată și cauzalitatea

# Packages and settings
packs <- c("ggplot2", "papaja")
success <- suppressWarnings(sapply(packs, require, character.only = TRUE))
install.packages(names(success)[!success])
sapply(names(success)[!success], require, character.only = TRUE)

theme_set(papaja::theme_apa())  # set APA style theme for plots

# Get data from package
if(!require(MultiKink)){install.packages("MultiKink")}; library(MultiKink)
data <- MultiKink::triceps
data <- subset(data, age <= 30)

# Plot the data
plot <- ggplot(data, aes(x = age, y = triceps)) +
  geom_point(alpha = .7, fill = "black", color = "grey50") + 
  xlab("Vârstă") +
  ylab("TSF")

# Piecewise regression
data_k1 <- data[data$age < 5, ]
data_k2 <- data[data$age >= 5 & data$age < 10, ]
data_k3 <- data[data$age >= 10 & data$age < 20, ]   
data_k4 <- data[data$age >= 20, ]
  
lm_p1 <- lm(triceps ~ age, data = data_k1)
lm_p2 <- lm(triceps ~ age, data = data_k2)
lm_p3 <- lm(triceps ~ age, data = data_k3)
lm_p4 <- lm(triceps ~ age, data = data_k4)

pred_p1 <- predict(lm_p1, se.fit = TRUE, data = data_k1)
pred_p2 <- predict(lm_p2, se.fit = TRUE, data = data_k2)
pred_p3 <- predict(lm_p3, se.fit = TRUE, data = data_k3)
pred_p4 <- predict(lm_p4, se.fit = TRUE, data = data_k4)

plot_pwr <- 
  plot + 
  geom_line(data = data_k1, aes(y = pred_p1$fit), size = 1, color = "red") +
  geom_ribbon(data = data_k1, aes(ymax = pred_p1$fit + pred_p1$se.fit, ymin = pred_p1$fit - pred_p1$se.fit),
              fill = "red", alpha = .1) + 
  geom_line(data = data_k2, aes(y = pred_p2$fit), size = 1, color = "red") +
  geom_ribbon(data = data_k2, aes(ymax = pred_p2$fit + pred_p2$se.fit, ymin = pred_p2$fit - pred_p2$se.fit),
              fill = "red", alpha = .1) +
  geom_line(data = data_k3, aes(y = pred_p3$fit), size = 1, color = "red") +
  geom_ribbon(data = data_k3, aes(ymax = pred_p3$fit + pred_p3$se.fit, ymin = pred_p3$fit - pred_p3$se.fit),
              fill = "red", alpha = .1) +
  geom_line(data = data_k4, aes(y = pred_p4$fit), size = 1, color = "red") +
  geom_ribbon(data = data_k4, aes(ymax = pred_p4$fit + pred_p4$se.fit, ymin = pred_p4$fit - pred_p4$se.fit),
              fill = "red", alpha = .1)  

suppressWarnings(print(plot_pwr))

# Packages and settings
packs <- c("ggplot2", "papaja", "osfr", "tidyverse", "scales", "lspline", "plotly")
success <- suppressWarnings(sapply(packs, require, character.only = TRUE))
install.packages(names(success)[!success])
sapply(names(success)[!success], require, character.only = TRUE)

theme_set(papaja::theme_apa())  # set APA style theme for plots

# Read data from folder -- here we will read directly from OSF repository
# csv_name <- "cognitive_dynamics_heino.csv"
# csv_path <- here::here("datasets", csv_name) 
# df <- read.csv(csv_path)   

# Read data from OSF repository
repo_url <- "https://osf.io/w9v28/"
csv_name <- "cognitive_dynamics_heino.csv"
repo <- osfr::osf_retrieve_node(repo_url); Sys.sleep(10)
repo_files <- osfr::osf_ls_files(repo)
# we download to a temporary file that gets deleted after
temp_dir <- tempdir()
osfr::osf_download(repo_files[repo_files$name == csv_name, ], path = temp_dir); Sys.sleep(10)
csv_path <- file.path(temp_dir, csv_name)
df <- read.csv(csv_path)
unlink(csv_path)   

# Transform
df_clean <- 
  df %>%
  # keep only days with record of date
  tidyr::drop_na(date) %>%
  # convert to ISO 8601 date format, format used is not consistent
  tidyr::separate(date, into = c("day", "month", "year"), sep = "\\.", remove = FALSE) %>%
  tidyr::separate(time_of_day, into = c("hour", "minute"), sep = "\\:", remove = FALSE) %>%
  # there are some typos, exclude them
  dplyr::mutate(across(c(day, month, year, hour, minute), as.numeric)) %>%
  dplyr::filter(year >= 2014, year <= 2017, month < 12, day < 31) %>%
  dplyr::mutate(
    date = lubridate::ymd(paste(year, month, day, sep = ' ')),
    date_posixct = as.POSIXct(date),
    hm = lubridate::hm(paste(hour, minute, sep = ' '))
  ) %>%
  dplyr::filter(date > lubridate::as_date("2014-07-01")) %>%
  dplyr::select(rowNumber, date, date_posixct, hm, stroop_incongruent_post_meditation) %>%
  tidyr::drop_na(stroop_incongruent_post_meditation) %>%
  dplyr::mutate(y = stroop_incongruent_post_meditation,
                moments = dplyr::row_number())

# Plot
plot <- 
  ggplot(df_clean, aes(x = date_posixct, y = y)) +
  geom_point() +
  scale_x_datetime(date_labels = "%Y-%m", breaks= scales::date_breaks("6 month")) +
  xlab("") +
  ylab("Timp răspuns Stroop inc. Post mindfulness") + 
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

# Packages and settings
packs <- c("tidyverse", "papaja", "latex2exp", "cowplot")
success <- suppressWarnings(sapply(packs, require, character.only = TRUE))
install.packages(names(success)[!success])
sapply(names(success)[!success], require, character.only = TRUE)

if(!require(ggbrace)){devtools::install_github("nicolash2/ggbrace")}; library(ggbrace)

theme_set(papaja::theme_apa())  # set APA style theme for plots


# Data
df <- tibble(
  y = c(c(3, 4, 6, 5, 4), c(8, 7, 8, 9, 9)),
  AB = c(rep("A", 5), rep("B", 5)), 
  t = 1:10
)

# Stats
df_A <- df[df$AB == "A", ]
df_B <- df[df$AB == "B", ]  
mean_A <- mean(df_A$y, na.rm = TRUE)
mean_B <- mean(df_B$y, na.rm = TRUE)
intercept_A <- summary(lm(y ~ t, data = df_A))$coefficients[1, 1]
slope_A <- summary(lm(y ~ t, data = df_A))$coefficients[2, 1]
intercept_B <- summary(lm(y ~ t, data = df_B))$coefficients[1, 1]
slope_B <- summary(lm(y ~ t, data = df_B))$coefficients[2, 1]

# Packages and settings
packs <- c("ggplot2", "papaja", "scda", "scan")
success <- suppressWarnings(sapply(packs, require, character.only = TRUE))
install.packages(names(success)[!success])
sapply(names(success)[!success], require, character.only = TRUE)

theme_set(papaja::theme_apa())  # set APA style theme for plots

# Get data
# Michael's data from `scda` package (can be downloaded from presentation Gihub)
data <- scda::Singh
data_Michael <- subset(data, id == "Michael")

# Plot
plot <- 
  ggplot() +
    geom_point(data = data_Michael, aes(x = time, y = score_physical)) +
    geom_rect(aes(xmin = 4.4, xmax = 4.6, ymin = 0, ymax = 3), 
              fill = "black", alpha = 0.2) +
    geom_line(data = data_Michael[data_Michael$phase == 0, ], aes(x = time, y = score_physical)) +
    geom_line(data = data_Michael[data_Michael$phase == 1, ], aes(x = time, y = score_physical)) +
    scale_x_continuous(breaks = 1:18, limits = c(1, 18)) +
    scale_y_continuous(breaks = 0:3, limits = c(0, 3)) +
    xlab("Timp") +
    ylab("Scor fizic") +
    coord_equal()